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Nature of complex singularities for the 2D Euler equation 
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A detailed study of complex-space singularities of the two-dimensional incompressible Euler equa- 
tion is performed in the short-time asymptotic regime when such singularities are very far from 
the real domain; this allows an exact recursive determination of arbitrarily many spatial Fourier 
coefficients. Using high-precision arithmetic we find that the Fourier coefficients of the stream 
function are given over more than two decades of wavenumbers by F(k) = C(9)k~ a e~ kS ^ e \ where 
k = k(cos0, sin9). The prefactor exponent a, typically between 5/2 and 8/3, is determined with an 
accuracy better than 0.01. It depends on the initial condition but not on 6. The vorticity diverges 
as s -13 , where a + fJ = 7/2 and s is the distance to the (complex) singular manifold. This new 
type of non-universal singularity is permitted by the strong reduction of nonlinearity (depletion) 
which is associated to incompressibility. Spectral calculations show that the scaling reported above 
persists well beyond the time of validity of the short-time asymptotics. A simple model in which 
the vorticity is treated as a passive scalar is shown analytically to have universal singularities with 
exponent a — 5/2. 



Und es wallet und siedet und brauset und zischt, 
Wie wenn Wasser mit Feuer sich mengt, 
Bis zum Himmel spritzet der dampfende Gischt, 
Und Flut auf Flut sich ohn' Ende drdngt... 

Fricdrich von Schiller, from Der Taucher [1] 



I. INTRODUCTION 

A quarter of a millenium has elapsed since Euler pub- 
lished for the first time what is now known as the Eu- 
ler equations of hydrodynamics [2] . There has not been 
much celebration but this may just reflect our embar- 
rassment at not having made enough progress. Actually, 
Leonhard Euler warned us. At the end of his 1755 paper 
he wrote: "However all that the Theory of fluids holds, 
is contained in the two equations above, so that in the 
pursuit of the research we are not lacking the principles 
of Mechanics, but solely the Analysis, which is not yet 
cultivated enough for this design: hence we see clearly, 
which discoveries are left for us to make in this Science, 
before we can attain a more perfect Theory of the motion 
of fluids" . (A paper in Latin Principia motus fluidorum, 
published a few years after the paper in French, contains 



1 In French: Cependant tout ce que la Theorie des fluides ren- 
ferme, est contenu dans les deux equations rapportees cy-dessus, 
de sorte que ce ne sont pas les principes de Mechanique qui nous 
manquent dans la poursuite de ces recherches, mais uniqucmcnt 
l'Analyse, qui n'est pas encore asses cultivee, pour ce dessein: et 
partant on voit clairement, quelles decouvertes nous restent en- 
core a faire dans cette Science, avant que nous puissions arriver 
a une Theorie plus parfaite du mouvement des fluides. 



the basic equations and was already presented under a 
different title to the Berlin Academy in 1752.) 

Euler considered both the compressible and incom- 
pressible cases. Here we are concerned only with the 
latter which is particularly difficult in view of the global 
nature of the incompressibility constraint. One of the 
most important open questions concerning the "analysis" 
of the Euler equations is the well-posedness: does ini- 
tially smooth 3D flow, which is known to remain smooth 
for short times, eventually "blow up" , that is become 
singular in a finite time (see, e.g. Refs. [3, 4])? In two 
dimensions it has been known since the thirties that flow 
in a bounded domain, initially sufficiently smooth, never 
blows up [5, 6]. It was also shown that if such a 2D flow is 
initially analytic it will stay so forever [7-9]. However, in 
the course of time, such flow can develop very fine scales 
and there is a large discrepancy between the analytic es- 
timation of how the smallest scale decreases in time (a 
double exponential) and what is found in numerical sim- 
ulations (a simple exponential; see, e.g., Ref. [10]). 

The likely cause of the discrepancy is depletion, the 
phenomenon by which high-Reynolds number or inviscid 
incompressible flow tends to organize itself into struc- 
tures having vastly reduced nonlinearities (see, e.g., 
Ref. [11]). Depletion, which is still very poorly under- 
stood, may hold the key for understanding why 3D high- 
Reynolds number flow seems never to blow up, at least 
in simulations. 2 In this paper we shall focus on the two- 



2 For the case of 3D inviscid Euler flow there is no truly conclusive 
evidence in favor of blow up [4, 17]. Furthermore, if the flow is 
initially analytic, any real singularity will have to be preceded 
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dimensional case. 

There are well-known 2D examples of depletion, such 
as flows which depend only on one Cartesian coordinate 
or on the radial polar coordinate. Such flows are how- 
ever steady and thus globally depleted, with no dynam- 
ics. In this paper we shall be interested in 2D flow with 
an initial stream function which is a real trigonometric 
polynomial in the space variables, of the sort already con- 
sidered in Rcfs. [4, 12]. These are the 2D counterparts 
of well-known 3D flows such as Taylor-Green and Kida- 
Pclz [13-17] which have been used for (so far inconclu- 
sive) investigation of finite-time blow up. Our 2D flows 
have generally non-trivial dynamics and display locally 
very strong depletion. 

Trigonometric polynomials are instances of entire func- 
tions, that is, functions which are analytic in the whole 
complex domain. The only singularities of such func- 
tions are at complex infinity. The solution of the Euler 
equations at times t > sufficiently small can then be ex- 
tended analytically into the complex domain [7-9] . There 
is strong numerical evidence in 2D and also in 3D that 
such flow does not stay entire and develops singularities 
at certain complex locations for any t > [4, 10, 12, 14]. 
Complex singularities are usually detected through the 
Fourier transforms of the solution: roughly, there is an 
exponential tail related to the distance of the nearest 
singularity from the real domain, accompanied by an al- 
gebraic prefactor related to the nature (also called type 
or structure). of the singularities. 

Little is known about the nature of complex singu- 
larities of the Euler equations. In Refs. [4] and [12] it 
is shown numerically for the 2D case with the initial 
stream function cosxi + cos 2^2 that the complex sin- 
gularities lie on a smooth manifold and that the vorticity 
becomes infinite when approaching the singular mani- 
fold; there is however considerable uncertainty as to the 
scaling law of this divergence. In Ref. [18] the motion of 
preexisting complex-space singularities is studied analyt- 
ically but their nature is kept quite arbitrary. In Ref. [19] 
traveling-wave solutions with a pure imaginary velocity 
are studied for 3D axisymmetrical flow with swirl; using 
an ultra-high precision 3 numerical method, the singular- 
ities in the complexified axial variable are mapped out 
as a function of the (real) radial variable and found to 
lie on a smooth curve; the nature of the singularities is 
obtained using a "sliding fit" method. In Ref. [20], for 
the vortex sheet problem with an initially analytic inter- 
face, the nature of complex singularities of the interface 
is obtained using an ultra-high precision method and a 
"pointwise fit" . The sliding fit and the pointwise fit are 
very closely related to the method we use in Section III A 
and we shall come back to this matter. Since the work 
of Krasny [21], it appears that ultra-high precision is a 



prerequisite for obtaining numerical information on the 
nature of singularities, particularly when they are in the 
complex domain. 

From a theoretical point of view, for many nonlinear 
equations of mathematical physics a very successful tool 
in studying the nature of singularities has been dominant 
balance and its refined versions such as Painleve analysis 
[22] . Dominant balance analysis typically gives universal 
singularities, that is singularities whose positions may de- 
pend on the initial conditions but their nature does not. 
The simplest instance is the ID viscous Burgers equa- 
tion whose complex-space singularities are simple poles, 
obtainable by balancing the nonlinear term against the 
viscous one. For the ci-dimensional incompressible Euler 
equations, attempts to use dominant balance fail because 
of the particular structure of the nonlinearity: if we as- 
sume that the solution becomes singular on a complex 
manifold of dimension d — 1 , the nonlinearity vanishes to 
leading order. This is just a consequence of the simplest 
form of depletion, the vanishing of nonlinearity for solu- 
tions which depend on a single spatial coordinate. The 
nature of singularities cannot be obtained by a dominant 
balance argument; actually, as we shall see, complex sin- 
gularities of the 2D Euler equation display a very unusual 
non-universality. 

In this paper we will mainly discuss the short-time 
asymptotic regime presented in Ref. [4] and extensively 
used in Ref. [12] which gives us the most accurate in- 
formation on the nature of complex singularities. 4 After 
briefly introducing it in Section II we will show that this 
regime can be reformulated as a "pseudo-hydrodynamic" 
Euler problem, in which all the action including the sin- 
gularities takes place in a plane extending in the pure 
imaginary directions, but our usual hydrodynamic intu- 
ition is still applicable. The short-time asymptotics al- 
lows us to obtain recurrence relations for spatial Fourier 
components involving only wavectors k = (k±, k 2 ) with 
fci > and k 2 > 0, a feature which is also present 
in the Moore approximation for vortex sheets [23] and 
its generalization to smooth flow [19]; as a consequence 
Fourier components can be calculated in ultra-high pre- 
cision without any truncation error. In Section III we 
present the numerical evidence for simple scaling laws 
associated to complex singularities and determine the 
nature of the singularities with high precision. Analyz- 
ing short-time asymptotics for different initial conditions, 
we find that the singularities are non-universal. In Sec- 
tion IV we describe the global and local geometry of the 
pseudo-hydrodynamic flow, including depletion of non- 
linearity which is especially strong near the singularities. 

Sections III and IV both involve a mixture of numeri- 
cal results and of theoretical arguments, some heuristic, 
some more rigorous. We must stress that at the moment 
we do understand various features of the solution, in par- 



by complex-space singularities [8, 9]. 
3 that is, higher than double precision. 



4 Henceforth, Ref. [12] will be cited as MBF. 
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ticular why the scaling exponent for singularities does not 
depend on the direction, but we failed so far to reproduce 
by theory the non-universal scaling exponents observed 
for the singularities. Nevertheless by moving to yet an- 
other level of toy-modeling (the equivalent for our prob- 
lem of considering the vorticity as a passive scalar in a 
prescribed velocity field), we can determine the nature of 
the corresponding complex singularities using dynamical 
systems tools (Section V). The nature of these "advec- 
tion" singularities is however universal and therefore does 
not reproduce an essential feature of the nonlinear Euler 
flow. Finally, conclusions, open problems and a tentative 
road map for future research on blow up are presented in 
Section VI. To make the present paper reasonably self- 
contained we shall occasionally re-derive results already 
found in Ref. [4] and MBF. 

II. SHORT-TIME ASYMPTOTICS AND 
PSEUDO-HYDRODYNAMICS 

We are interested in the short-time asymptotics for 
the 2D Euler equation, written in terms of the stream 
function 

d t \7 2 ^(x, t) - J(*, V 2 *) = , (1) 

where x = (x 1 ,x 2 ) and J(f,g) = difd 2 g - d\gd 2 f '. 
The initial condition ^?q{x) = ^>(x,0) is a real 2n- 
periodic trigonometric polynomial of the form *S?o( x ) = 
Y /k F ( -^(k)e ik - x , where the sum has only a finite number 
of terms. Here k — (ki,k 2 ), where fci and fc 2 are signed 
integers. The short-time asymptotics is simplest when 
the initial condition has only two orthogonal Fourier 
modes, as in Refs. [4, 12] where the assumed initial con- 
dition is 

^o(x) = cosxi +cos2a;2 • (2) 

In what follows we shall call this initial condition Stan- 
dard Orthogonal Case (SOC). One of our present goals 
is to investigate to what extent complex singularities are 
or are not universal, we are thus naturally led to con- 
sidering more general cases, having, e.g., more than two 
modes in the initial conditions. In Appendix A it will 
be shown that the short-time asymptotic regime for the 
multimode case can be reduced to a set of two-mode ini- 
tial conditions. We may thus without loss of generality 
limit ourselves to two-mode initial conditions of the form 

y (x) = h lC ip a: + h 2 c iq x + c.c. . (3) 

Here c.c. stands for "complex conjugate", p — (pi,p 2 ) 
and q = (qi, q 2 ) are two vectors with signed integer com- 
ponents. Furthermore, we assume that p and q are not 
parallel and do not have the same modulus since other- 
wise the two-mode initial condition is a time- independent 
solution of the Euler equation. By performing if needed 
a suitable translation, we can then assume that hi and 



h 2 are positive. Finally, since our goal here is primarily 
to demonstrate non-universality of the nature of the sin- 
gularities with respect to the initial conditions, we shall 
not strive for the greatest generality and limit ourselves 
to basic modes with non-negative components such that 
Piq2 - qip2 > 0. 

Eq. (1) has a solution in the form of a Taylor series in 
the time variable 

*(x,t) = ^* n (x)t n , (4) 

n>0 

where is the initial condition and the *„(a;)'s for n > 
1 are easily shown to satisfy the recursion relations: 

V 2 *„ +1 = — ^- J2 J(y m y 2 V P ). (5) 

va+p—n 

For the two- mode initial condition (2) all the ^ n {x) are 
trigonometric polynomials that can be continued analyt- 
ically to complex locations z = x + iy. Since the initial 
condition has its singularities at infinity, we expect, by 
continuity, that at short times the singularities will have 
large imaginary parts \yi\ and \y 2 \. Let us now suppose 
that 2/i — > +oo and y 2 — > +oo in such a way that their 
ratio y 2 /yi stays finite but arbitrary. Obviously, the four 
vectors (pi,p 2 ), (-pi,-p2), (qi,q2) and (-qi,-q 2 ) di- 
vide the fe-space into four angular sectors so that, e.g. in 
the first angular sector P2/P1 < y 2 /yi < q2/qi- Then, for 
z such that y lies in the first angular sector, to leading 
order any additional factor t in the expansion (4) is ac- 
companied by cither a factor e _Ip z or a factor e~ iqz , 
thus giving amplitude factors te py and te qy , respec- 
tively. When t — > and \y\ — > 00 these factors remain 
finite, provided p ■ y and q • y are shifted by Int. This 
suggests that the short-time asymptotics is obtained by 
the similarity ansatz 

*(z,t) = (l/t)F(z), (6) 
z = (zi, z 2 ) = (zi -HAilnt, z 2 + iA 2 lni) , (7) 

where Ai and A2 are determined by 



n\ — ciuu /\ 2 — • 1UI 

qip2 - q2pi qip2-q2pi 
Substitution in (1) gives the similarity equation 

V 2 (-l + i XA, +i\ 2 d Z2 )F = J(F, V 2 F) , (9) 

where the tilde means that the partial derivatives are 
taken with respect to the new variables. The initial con- 
dition (2) becomes an asymptotic boundary condition 

F(z) ~ hie~ lp z + h 2 e~ ict ' z , y\ — > — 00, y 2 — > -00 . 

(10) 

In (9) the second and third terms on the l.h.s. can be 
viewed as stemming from the advection by a pure imag- 
inary constant "drift velocity". This is because we are 
following the singularities coming "down" from complex 
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infinity. It is important to observe that (9) is an exact 
consequence of the Euler equation. The only place where 
an approximation is made is in the boundary condition 
(10) where harmonics containing e.g. c +lp z and c +iq z 
are discarded because such terms are exponentially sub- 
dominant at short times. 

In what follows we shall generally limit ourselves to 
SOC, giving occasionally an indication of what is valid 
for more general two-mode cases. The general case can 
easily be handled but we wish to avoid burdening the 
reader with unnecessarily complicated statements and 
equations. 

The function F(z), which is 27r-pcriodic in X\ and 7r- 
periodic in x 2 , is analytic in the product of the half-spaces 
yi < and y 2 < and thus its spatial Fourier series has 
only harmonics of the form e _1 ( felZl+fe2Z2 ) with k\ > 
and fc 2 > 0. Its Fourier series is here written as 5 

oo oo 

F{z 1 ,z 2 )= ^(-l) fel ^(fci,fc 2 )e- ifclil e- ifc2i2 . 
k 1= a k 2 =a 

(11) 

The reason for the presence of the factor (— l) kl will 
become clear shortly. The Fourier coefficients F(k) = 
F(ki,k 2 ) can be calculated recursively from the relation 
given in MBF which follows from (9) 

P ^ = - k 1 + k 2 ,2-l W X (12) 
fci fe 2 

£(pAk)|k-p| 2 F(p 1 ,p 2 )F(fc 1 -p 1 ,fc 2 -p 2 ) . 

Pl=0p 2 =0 

Here, p A k = p\k 2 — p 2 k\. Because there are no Fourier 
harmonics with negative fci or fc 2 , the convolutions in (12) 
only involve positive arguments. This feature, which al- 
lows truncation-free determination of Fourier coefficients, 
is also present in the Moore approximation for the vortex 
sheet problem and in its generalization to axisymmetri- 
cal flow [19, 23]. The initialization of the recursion rela- 
tions requires the knowledge of the coefficients along the 
"edges" , that is the half-lines fci = and fc 2 = 0. In the 
present case F(0,2) = 1/2 and F(1,0) = -1/2 while all 
the other edge harmonics are zero. 6 It has been shown in 
MBF that, with the choice made above in (11), the co- 
efficient F(1,0) is the only one that is negative. All the 
other ones are non-negative. This result has so far only 
been established by (very solid) numerical computations 
and holds for all the two-mode initial conditions studied. 
As we shall see, this has important consequences for the 
geometry of singularities. 



We shall now show that (9) can be reformulated as 
the steady solution of a pseudo-hydrodynamic problem 
in a suitable imaginary plane. Since we are working with 
analytic functions, we can replace the complex partial 
derivatives d Zl and d Z2 by — id Vl and — id y2 , holding the 
^-coordinates fixed. In terms of such y-derivatives (9) 
becomes an equation with real coefficients. If we fur- 
thermore choose xi and x 2 such that the boundary con- 
dition (10) becomes real then the solution "above such 
points" F is also real. This happens for x\ = 0, n and 
for x 2 = 0, 7r/2, 7r, 37r/2. The positivity of all but one 
of the Fourier coefficients defined in (11) with the factor 
(— 1) , amounts to stating that, after moving the origin 
to (ir, 0), all but one of the usual Fourier coefficients of F 
are positive. As we shall see in Section IV, this gives us 
the possibility of analyzing the (short-time) complex sin- 
gularities by focusing solely on the y-plane above (ir, 0). 
This point turns out also to be a center of symmetry for 
the Euler flow with the initial condition (2), but it is not 
clear whether this matters. 7 Henceforth we shall consider 
the j/-plane above (ir, 0). 

We define a pseudo-stream function in terms of the 
y-coordinates (from now on we drop the tilde on the y 
variables for notational simplicity) 

i>{y) = i>{yi,V2) = ^yi -y 2 + F(tt + iy u iy 2 ) ,(13) 

= \ yi -y 2 +J2Hk)c k v , (14) 
k 

where the two linear terms on the r.h.s. have been in- 
troduced to avoid having an additional advection term. 
Note that because of these terms, -0 is not the continua- 
tion to complex coordinates of a function periodic in x\ 
and x 2 . 

It is now elementary to check that (9)-(10) are equiva- 
lent to taking the steady-state (r-independent) solution 
of the pseudo-hydrodynamic equation 

<9 T V 2 V - J(V>, V 2 ^) = - V V , (15) 

with the asymptotic boundary condition (for yi, y 2 — > 

— oo) 

HvuV2) - \vi +V2^ Af 1 + \? m ■ (16) 

Here, in order to bring out familiar hydrodynamic nota- 
tion, we have introduced a pseudo-time variable t. 8 Wc 
are using V = (d\, d 2 ) for V y and the Jacobian J has 
its usual definition in terms of y-derivatives. We now in- 
troduce a pseudo-velocity and a pseudo-vorticity by the 



5 In MBF fci and &2 were denned with the opposite sign. 

6 In the general case of two basic modes p and q, the main change 
with respect to SOC is the replacement in (12) of the denomina- 
tor fci +&2/2 — 1 by Aifci + A2&2 — 1 where Ai and A2 are defined 
in (8). 



Note that streamlines have a hyperbolic structure near (it, 0), 
but an elliptic structure near (0, 0) which is also a center of 
symmetry. 

If we allow the function F and thus tfi to also depend on t and 
set r = ln(l/t), we obtain precisely (15). 
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usual definitions: 9 

v = (ui, v 2 ) = (d 2 ip, -<9iV0 (17) 

= -(1, 1/2) + ]T (k 2 , -fci) F(k) e k y , (18) 
k 

uj = -V 2 V, (19) 

= H - fc2 ^( fc ) cfe y . ( 2 °) 
fe 

in terms of which (15) reads 

[d T oj] + v ■ Vlu + lo = , (21) 
with the boundary conditions (for yi, y 2 — > — oo) 

„ ~ le" - 2e 2 «* . (22) 

For other initial conditions, only the boundary condi- 
tion (22) must be modified. The T-derivative term has 
been put within square brackets since we are only inter- 
ested in the steady-state solution. Note that the pseudo- 
hydrodynamic formulation in the y-plane is that of a 
quasi-two-dimensional flow in a 3D container with bot- 
tom friction producing a Rayleigh drag. In this formu- 
lation r — > +oo as we approach the initial instant. An 
alternative interpretation is to define r as In t, to avoid 
reversing the course of time, and then to change the signs 
of v and of to and replace the Rayleigh drag by an insta- 
bility. 

In the pseudo-hydrodynamic formulation it is now ob- 
vious that the problem is invariant under an arbitrary 
translation h = (hi, h 2 ) in y-space. By (14), such a 
translation amounts to a factor e k h on the Fourier co- 
efficients F(k). It follows, as noted in MBF, that the 
set of initial conditions ^fo( x ) = ° hl cos mi + e 2/l2 cos 2x 2 
is equivalent to SOC as long as h is within the analyt- 
icity domain. Similarly, a translation in fc-space with 
integer components (ni, n 2 ) is equivalent to multiplying 
F(n + i 3/i, iy 2 ) by the exponential factor e niVl+ri2V2 in 
y-space. The exponential being an entire function, this 
changes neither the positions nor the nature of the sin- 
gularities at finite distance. 

III. NUMERICAL INVESTIGATION OF 
SCALING LAWS IN FOURIER SPACE 

We shall show in this section that the solution of the 
Euler equation in the short-time asymptotic regime de- 
fined in the previous section, has remarkably clean scal- 
ing properties in Fourier space. By this we mean that 
the wavenumber dependence of the Fourier coefficients 



The true velocity is actually pure imaginary in the i/-plane and 
the true vorticity is —ui. 



is represented as a decreasing exponential multiplied by 
an algebraic prefactor whose exponent can be measured 
very accurately. Such a functional form is not surprising. 
In fact the exponential is the signature of the location 
of a singularity while the prefactor encodes the nature of 
the singularity. For one-dimensional analytical functions 
with isolated singularities in the complex space this is 
well known: a singularity at z + of the form (z — z+) p has 
a signature in the modulus of the Fourier transform at 
high wavenumbers k of the form C|/c|~ p ~ 1 c~ ,5 l fe , where 
S is the the distance of to the real axis (see, e.g., 
Ref. [24]). Such asymptotic results have been extended 
in the nineties to the Fourier transforms of periodic an- 
alytical functions of several complex variables when the 
wavevector k tends to infinity with a fixed rational slope 
tan 9 = k 2 /ki — p/q, where p and q are relative prime 
integers [25-27]. 

When the Fourier coefficients are obtained numerically, 
there is a maximum wavenumber /c max . Unless it is taken 
very large, there will be very few points on the line of 
slope p/q as soon as q is not a very small integer. But 
a large value of fc max entails extremely small Fourier co- 
efficients because of the exponential decrease with the 
wavenumber. Thus, as stressed in MBF, very high preci- 
sion may be needed to avoid swamping by the rounding 
errors. Truncation errors are not an issue in the short- 
time asymptotic regime since the Fourier coefficients can 
be calculated from (12) with arbitrary accuracy. 

The data obtained for the SOC initial condition in 
MBF had wavenumbers k= |fc| up to 1000 or 2000, de- 
pending on the direction and were calculated with 35- 
digit accuracy. 10 Most of the results presented here are 
based on the 35-digit calculation. Additional calcula- 
tions are also presented here with various initial condi- 
tions, with up to 100-digit precision and wavenumbers 
which can reach 4000 in particular directions. We note 
that the MPFUN90 package for high-precision calcula- 
tion used in MBF, here and in Ref. [19] makes use of fast 
Fourier transform techniques. Thus the CPU time per 
multiplication, as a function of the number of digits N, 
is proportional to N log TV [28] . 

We now show that it is quite easy in principle to ob- 
serve scaling by analyzing the behavior of the Fourier 
coefficients in directions of rational slope. Figs. 1 and 2 
give two examples of the analysis of Fourier coefficients 
along straight lines through the origin 11 in a direction of 
rational slope, using the data from MBF for the Fourier 
coefficients of the stream function with SOC initial con- 



10 In MBF it was stated that, when using only double-precision (15- 
digit) accuracy, unacceptably large errors are obtained beyond 
wavenumber 800. Actually, as pointed out by P. Zimmcrmann 
(private communication), the double-precision calculation can be 
modified in such a way that, up to wavenumber 1000, the relative 
error on Fourier modes does not exceed 10 — 5 . 

11 All the lattice lines of a given rational slope have the same high-fc 
asymptotics, due to the observation made at the end of Section II. 
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FIG. 1: Fourier coefficients of the stream function F along 
two lines of different slopes as a function of k = \k\ in lin-log 
coordinates. 
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FIG. 2: Same as in Fig. 1 after division by exp(— 5k) (com- 
pensated Fourier coefficients) in log-log coordinates. Most of 
the points are in the asymptotic power-law regime, at least 
visually. 



ditions. The first case has k^/hi — 2, the direction with 
the largest number of grid points having non-vanishing 
Fourier coefficients. The second case has A^/fci = 18/11, 
the direction with the slowest decrease of the Fourier co- 
efficients. Fig. 1, which shows the Fourier coefficients 
in lin-log coordinates,reveals an exponential tail oc e~ Sk ; 
a least square fit gives 5 = 0.021 for the first case and 
5 = 0.0065 for the second case. 12 In Fig. 2 we show the 
"compensated" Fourier coefficients obtained by dividing 
by the exponential term; the result is then represented 
in log-log coordinates in order to look for an algebraic 



12 Why the minimum value is so small is a matter we shall come 
back to in Section V. 



prefactor oc k~ a . The quality of the scaling obtained is 
impressive: over most of the range we cannot on a log-log 
plot visually distinguish the prefactor from a power law 
with exponent a = 8/3. As we shall see, the exponent 
does not depend on the direction chosen. 



A. Technique for capturing algebraic prefactors 

Determining the scaling properties as done above by 
use of least square fits, compensating exponentials and 
log-log plots is not optimally adapted for delicate issues 
such as studying the dependence of the prefactor expo- 
nent on the direction of the wavector or on the initial 
conditions. As pointed out by Shelley [20], it is better to 
remove some of the subjective biases present in a least 
square fit (such as choosing the range in k). We shall 
make use of his method of point-wise fit (also used in 
Ref. [19], where it is called a sliding fit), followed by an 
extrapolation step as now explained. 

In fc-space, a direction of rational positive slope is char- 
acterized by fo/fti = tan0 = q/p (where the positive in- 
tegers p and q are taken to be relative primes). All the 
k vectors on the line of slope p/q through the origin are 
thus of the form k = nko, where feo = (p, q) and n is a 
positive integer. What we have seen at the beginning of 
Section III suggests that for a given direction of rational 
slope tan#, the Fourier coefficients of the stream function 
can be represented, for sufficiently large k, as 



F ~ C(9)k 



- a (9)„-5(0)k 



(23) 



Henceforth a, C and S will be referred to as the prefactor 
exponent, the constant and the decrement, respectively. 
When there is no ambiguity, the 6*-dependence will be 
omitted. Following Ref. [20], let us assume for a moment 
that (23) holds exactly and let us set F n (k ) = F(nk ). 
It then follows that if we know F n (k ) for any three con- 
secutive values, say n — 1, n and n+1, we can determine 
C, a and 5 by 



In 



a = 



f 1 n-l(fc )f 1 » + i(fco) 



In 



(n - l)(n + 1) 



S = 



1 




(24) 



(25) 



InC = lnF„(fc ) + oln[(n)|fc |] + n\k \S . (26) 

The expression (24) for a follows immediately by notic- 
ing that in the combination F„_iF„ + i/F 2 the constant 
C and the exponential factor both drop out. The other 
two expressions are readily established by taking the log- 
arithm of (23). 
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FIG. 3: Local prefactor exponent ai oc (A;) vs wavenumber for 
two values of the slope. 
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FIG. 4: Discrepancy between 15- and 35-digit calculation of 
the local prefactor exponent a\ oc (k) along fe/fci = 5/3, as an 
estimate of the absolute error on ai oc - 



Of course, we have no reason to expect that (23) holds 
exactly for arbitrary wavenumbcrs. At best it will hold 
asymptotically at large wavenumbcrs. Nevertheless we 
can use (25)-(26) to calculate a local prefactor exponent 
ctioc(k), which depends on the wavenumber. Here we 
have chosen to use as arguments of the local quantities 
the wavenumber k = n|feo|- 

The typical behavior of a\ oc (k) is shown in Fig. 3 for 
two directions. For large values of k the curves grow to an 
asymptotic value close to 8/3. Globally, a\ oc (k) is found 
to be non-monotonic when 9 < 9* with tan^* close to 3 
(but not very sharply defined) and monotonic above 9*. 

To estimate the asymptotic value we must extrap- 
olate the data beyond the largest available wavenumber 
at which they are known with acceptable accuracy. Since 
the only cause of error in a\ oc (k) are rounding errors, we 
can measure such errors by comparing runs having dif- 
ferent levels of precision. Fig. 4 shows the discrepancy 



(absolute error) of a\ oc (k) obtained with 15 and 35 digit 
precision. The error is seen to grow with the wavenum- 
ber in an approximately exponential fashion, the highest 
value being about 10~ 8 around wavenumber 1000. We 
shall see that the error involved in the extrapolation may 
be much larger than 10 -8 . 

One well-known difficulty with extrapolation is that 
the problem may not be well-posed unless one has addi- 
tional information on the functional form of the conver- 
gence to zero of the remainder — a\ oc (k). In Ref. [20], 
which deals with the shape of a vortex sheet continued 
analytically to complex parameters, it is assumed that 
branch singularities of unknown exponent are present and 
that the high-fc behavior of the one-dimensional Fourier 
transform can be obtained from Laplace's method to 
leading and first sublcading orders; the inclusion of the 
first subleading correction allows a much improved de- 
termination of the exponent. This extrapolation pro- 
cedure is equivalent to assuming that the remainder 
®oo — ct\ oc (k) goes to zero as 1/k. For our problem, 
unfortunately no simple functional form of the remain- 
der, such as algebraic, exponential or inverse logarithmic 
decrease, gives a satisfactory fit. An efficient extrapo- 
lation method for a wide range of functional behaviors 
of the remainder is the epsilon algorithm of Wynn [29] , 
related to the Shanks transform method [30]. It is an 
algorithm for acceleration of convergence of a sequence 
S = (s(°), s< 2 ), sW) <E C, and it comprises 

the following initialization and iterative phases: 
Initialization: For n = 0, 1, 2, . . . 



e W = (artificially), = a(n) 

Iteration: For n = 0, 1, 2, . . . 



-(«) 
'l+i 



-(n+l) 



+ 



_(n+l) 



-(«) 



(27) 



(28) 



After a few iterations of the algorithm, applied to 35- 

(n) 

digit SOC data, the s] s with even I become almost 
constant and give an estimate of the extrapolated expo- 
nents (see Fig. 5). The epsilon-algorithm extrapolated 
exponents will be used when discussing results (unless 
otherwise stated). We have also used the recently intro- 
duced asymptotic interpolation method of van der Ho- 
even [31] which strips off successively leading and sub- 
leading terms by suitable transformations before doing 
the interpolation. This method works impressively for 
the passive scalar model discussed in Section V for which 
both leading and subleading terms in the high-fc expan- 
sion can be determined from numerical data. In the non- 
linear case, the asymptotic interpolation method gives 
exponents consistent with those determined by the ep- 
silon algorithm with a relative error of about 10~ 3 ; we 
have so far not been able to determine numerically the 
functional form of sublcading corrections. As we shall sec 
in Section IV C, theory tells us that a should not depend 
on the angle 9. We suspect that ^-dependent subleading 
corrections account for the slight apparent variation of a 
with 9, reported in Section IIIB. 
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FIG. 5: Local prefactor exponent a\ oc (n) for the nth point FIG. 7: Angular dependence of the decrement 6(6) for SOC. 
along the line fe/fci = 5/3 which has (ki, fe) = (5n, 3n); 
it is shown together with its second and fourth order epsilon- 
algorithm extrapolated values. Inset: enlargement for n > 50. 
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B. Results for SOC 



For the SOC, whose initial stream function is cosxi + 
cos 2x2, we now use the method described in Section III A 
to calculate the prefactor exponent a(6), the decrement 
6(6) and the constant C(6). 

Figs. 6, 7 and 8 show the angular variation of a, 6 
and C, respectively, excluding near-edge ranges where 6 
is close to or tt/2 which deserve separate discussion (see 
Section HID). 
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FIG. 8: Angular dependence of the constant C(6) for SOC. 
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FIG. 6: Angular dependence of the prefactor exponent a(6) 
for SOC (extrapolated by the epsilon algorithm) . Below 6 = 
0.2-7T (long dashed line) the extrapolation cannot be trusted. 



The most striking result is the very weak angular de- 
pendence of the prefactor exponent, which over the range 
0.2tt < 6 < 0.45?r is given by a = 2.66 ± 0.01, con- 
sistent with the theory which predicts independence on 
9 (Section IV C). This immediately leads to asking if 
asoc = 8/3. The short answer is: we do not know. We 
shall come back to this at length. 

The angular dependence of 6 has already been reported 
in MBF where it was measured by decomposing the set of 
directions into small angular sectors. 13 We find that 6(8) 
achieves a minimum value 6* « 0.0065 at 8* 0.3247T 
and it becomes large near the edges. In MBF it was 



In MBF 8 was varying in the third quadrant; here, because of the 
aforementioned change of notation 8 varies in the first quadrant. 
Furthermore, the method used in MBF was less accurate than 
the present one and there are thus small discrepancies in the 
values reported. 
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reported that the shell-summed amplitude of F(k) 

A(k)= £ |F(fe)|, (29) 
fc<|fe|<fc+i 

a kind of discrete angle average, behaves as C'k~ 216 e~ s * k 
for large k. This is consistent with the present result. 
Indeed for large k, we can evaluate the shell sums (29) by 
integrating over kdd using (23) and steepest descent near 
6+. This changes the prefactor from k~ a to fc-^+i- 1 / 2 ~ 
fc- 2 - 16 . 

Finally, C(9) is quite flat in the interval 0.1 < 9 < 0.4. 



C. Non-universality of the scaling exponent 

Having established the angle-independence of the pref- 
actor exponent, we now investigate its dependence on the 
initial condition. What happens when we change from 
SOC (given by (2)) to another initial condition? Since 
35-digit computations take up to one month of CPU, we 
generally used 15-digit accuracy but there is one impor- 
tant exception (see below). At first we changed SOC to 

^o(x) = cosxi + cos3x 2 , (30) 

for which the basic modes in the short-time asymptotics 
are (1, 0) and (0, 3) between which there is the same 
90-degree angle as for SOC. The prefactor exponent was 
again indistinguishably close to 8/3. For a while this led 
us to conjecturing the universality of the 8/3 exponent. 
Well . . . until we tried 

^o(x) — cos(xi + x 2 ) + cos 2x2 , (31) 

whose basic modes are (1, 1) and (0, 2), forming an angle 
of 45 degrees. This gave us an exponent a w 2.54. The 
same exponent was obtained with 

'i'o(x) = cos(xi + x 2 ) + cospx 2 , (32) 

with p = 1, 3, 4, whose basic modes are different but 
also form an angle of 45 degrees. We also did some ex- 
ploration of the direction dependence of a and, just as for 
SOC, did not find any. As we shall see in Section IV C, 
independence on the direction can be shown to hold. 

All this was pointing towards non-universality of the 
prefactor exponent, that is dependence on the initial con- 
dition or at least on the angle between the basic modes. 
To ascertain the non-universality we performed a 100- 
digit computation for (31) with /c max = 1000. Fig. 9 
gives the epsilon-algorithm extrapolated values of the lo- 
cal prefactor exponent for this calculation as a function 

of e. 
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FIG. 9: Angular dependence of the prefactor exponent a(6) 
for the "45-degree" initial condition ^o{x) = cos2a;i + 
cos(a;i + X2) 

Except near the edges the exponent stays very close 
to 2. 54. 14 The discrepancy between 2.54 and 2.66 vastly 
exceeds the estimated error on the prefactor exponent, 
as discussed in Section III A. Finally, we report that for 
all cases discussed in this section on non-universality, the 
positivity of all the Fourier coefficients except one, holds 
just as for SOC. 



D. Intermediate asymptotics near the edges 

In this Section we discuss only SOC, but the theoret- 
ical results presented are easily generalized. We have 
seen that on any line of strictly positive and finite ratio- 
nal slope the Fourier coefficients decrease exponentially 
at high k (up to algebraic prefactors). This is not true 
for lines of vanishing and infinite slope. We can explicitly 
calculate from the recursion relation (12) all the coeffi- 
cients having either k 2 = 2 or k\ = 1. Indeed along 
such "edge lines" the recursion relations take the form of 
first-order linear homogeneous finite difference equations 

1 h 2 — Oh, -A- A 

XL*.)- 5*^1^ *!.*»-»>■ (34) 

At large orders, essentially each coefficient on a horizon- 
tal or vertical edge line is obtained by dividing by k\ or 
fc 2 /2 the adjacent lower order coefficient. Thus they are 
decreasing roughly as 1/fci! or l/(fc 2 /2)!. More precisely, 
using standard asymptotic methods for difference equa- 



The anomalously low value around 8/n = 0.42 is caused by a 
large denominator in the corresponding slope (79/21) which docs 
not permit a reliable determination of a. 
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tions [32] , it is easily shown that for integer m — > oo 

F(m, 2) ~ F(l, 2m) - m^V 71 ™-" 1 , (35) 

which decreases faster than exponentially. 

If we now consider a "near edge" direction with 9 close 
to or to tt/2 we expect that the edge behavior will man- 
ifest itself as intermediate asymptotics making it hard 
to obtain clean scaling for the prefactor. We can how- 
ever easily predict the ^-dependence of the decrement 5 
by the following argument. When 9 is small, the line 
through the origin of slope tan 9 w 9 will intersect the 
edge &2 = 2 at fci w 2/9. At this point, by (35), the log- 
arithm of the Fourier amplitude is given to leading order 
by —(2/9) ln(2/0). Assuming that, on the line of slope 9, 
this point is within the region of exponential fall-off with 
decrement 5(9), we obtain 

-(2/0)ln(2/0) « -(2/9)5(9) , (36) 

which gives (for 9 — ► 0) 

In . (37) 

Near the other edge, we obtain by a similar argument 
(for 9 -> tt/2) 

^4 h (^)- (38) 

We turn now to numerical study of the near edge be- 
havior of Fourier coefficients. So far we have determined 
such coefficients in regions having comparable extensions 
in the k\ and k 2 directions. The structure of the recursion 
relation allows us however to determine the coefficients 
in rectangular domains having a very small or very large 
aspect ratio. We have seen that the local prefactor expo- 
nent behaves non-monotonically with the the wavenum- 
ber when 9 is below a critical value. 15 Consistently, we 
have found that for small 6>'s a more complex behavior is 
observed than for 0's close to ir/2. We have thus studied 
the former in more detail. Because of the slow conver- 
gence to asymptotics we need wavenumbers much larger 
than in MBF, so we used rectangular domains of size 
4000 x 480 near 9 = and of size 200 x 4000 near 9 = tt/2. 
Fig. 10 shows the variation with the wavenumber of the 
local prefactor exponent a\ oc (k) for various small 6's. It 
is seen that when 9 decreases, the wavenumber at which 
ttioc(fc) achieves its minimum increases and thus the ex- 
trapolation of a becomes more difficult. The situation is 
much more favorable for the determination of the decre- 
ment, because it is (logarithmically) large. 



This may be related to the fact that the global structure seen in 
Fig. 14 is far from being symmetrical in j/i and j/2- 
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FIG. 10: Wavenumber dependence of local prefactor exponent 
aioc for various small 6. Inset: da ioc (k)/dk. 






3(9) + 
ln(2/6)-l _ 
















."Ssfc 








0.1 

9 



FIG. 11: Angle dependence of the decrement S for small 9 
in lin-log coordinates (crosses). The continuous line is the 
theoretical prediction. 

Fig. 11 shows the measured decrement together with a 
theoretical prediction 5(9) = ln(2/0) — 1 which includes 
a sublcading correction to the leading-order prediction 
(37), obtained by a partially heuristic procedure. Near 
9 = tt/2 the decrement has also logarithmic scaling (not 
shown) , consistent with the leading-order prediction (38) 
but not very clean. As to the constant C(9), we found 
that it becomes large near the edges. For 9 — > the 
behavior is roughly C(9) cx 1/9 but there are substantial 
uncertainties because the constant C is quite sensitive to 
small errors made on 5 and a. 



E. Beyond short times 

In MBF it was shown, for the SOC initial condition, 
that deviations from short-time asymptotics become im- 
portant around t = 0.1. More precisely, deviations from 
the law 5(t) cx ln(l/t) become visible (see Fig. 2 of MBF). 
We now investigate numerically the issue of persistence 
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FIG. 13: Absolute value of the Fourier coefficients of the 
stream function at t = 0.8 along the rational direction 
kz/ki — 1 in lin-log coordinates. A least square fit (con- 
tinuous line) gives cfc _2,68 e _0 ' 429fe . The inset shows the same 
data after division by e -0 ' 429fe in log-log coordinates. 



FIG. 12: Contours of the absolute value of the Fourier coef- 
ficients (logarithmic scale) of the stream function at t — 0.8 
by full spectral simulation for SOC. 

of the k~ 2m law for SOC beyond the time of validity of 
short-time asymptotics. For this we must use a full spec- 
tral simulation with time-marching as in Refs. [4, 12]. A 
priori there is no need to use a resolution in excess of 
1024 2 since we shall see that the k~ 2m law deteriorates 
significantly after t = 1. At that time, the decrement 
S w 0.4, which implies that the flow is extremely well 
resolved with 1024 2 modes. Fig. 12 shows the behavior 
of the absolute value 16 of the Fourier coefficients of the 
stream function at t = 0.8 in the (k\, fc 2 )-plane (because 
of the Hermitian symmetry we are not showing negative 
ki). It is seen that there is a direction of slowest de- 
crease which has k 2 /k\ ~ 1- At short times the slowest 
decrease had k 2 jk\ k> 18/11 but this direction changes in 
the course of time. Fig. 13 shows the usual exponential 
decrease with an algebraic prefactor for the Fourier am- 
plitude in the direction of slope unity at t — 0.8. Beyond 
wavenumber 85, rounding errors take over (the calcula- 
tion has 15-digit precision). The same procedure, applied 
at much later times, e.g. at t = 1.9, still shows some kind 
of exponential tail but the data are far too wiggly to per- 
mit the extraction of a reliable power-law prefactor. We 
have also repeated the analysis beyond short times for the 
flow with initial condition cos(xi + x 2 ) + cos2x 2 , where 
the basic modes make an angle of 45 degrees. At time 
t = 1.4 the prefactor exponent is around 2.58, quite close 



Because of the symmetry of SOC, the Fourier coefficients are 
real. 



to the the value 2.54 reported at short times. 

Let us now discuss some of the limitations involved 
in the search for prefactor scaling beyond the short-time 
asymptotics. We begin with practical limitations. With 
1024 2 modes the direction of rational slope k 2 /k\ = 1 
has about 30 points before encountering 15-digit round- 
ing level. Other directions have typically only ten points 
and this makes precise determination of the decrement 6 
and the prefactor exponent a impossible. We thus can- 
not comment on any possible angular dependence of a. 
Calculations with higher resolution require higher preci- 
sion in order to lower the rounding noise level and this 
in turn requires enormous computer resources by a time- 
marching full spectral method if we demand that tempo- 
ral truncation error be at rounding level. 

There is a more fundamental issue regarding the va- 
lidity of the short-time asymptotic regime. For SOC, 
this regime breaks down around t = 0.1, as far as the 
temporal behavior of 8{t) is concerned. Actually the 
short-time approximation is strongly non-uniform with 
respect to the wavenumber: high wavenumbers show dis- 
crepancies at much earlier times than 0.1. For example 
we know that in the short-time regime, all the Fourier 
coefficients except one are non-negative, but as early as 
t = 10~ 3 , a 90-digit calculation by time-marching shows 
that Fourier coefficients start oscillating in sign beyond 
wavenumber forty. 17 . By t — 0.8 such oscillations arc 
found in the 15-digit calculation whenever k 2 jk\ > 2, 



It matters how precisely we let t — ► and k — + oo. For a fixed 
value of t, however small, the high-fc regime discussed in most of 
this paper may be just an intermediate asymptotic regime. It is 
conceivable that the non-universality found here is confined to 
this particular asymptotic regime. 
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irrespective of wavenumber. In the presence of such os- 
cillations, the functional form we have used in the short- 
time asymptotics oc k~ a e~ Sk is clearly invalid. What is 
happening has a geometric interpretation which is more 
readily understood after reading the first page of Sec- 
tion IV. In the short-time regime the positivity of the 
Fourier coefficients implies that the singular manifold is 
in the y-plane. Note, however, that in this regime we are 
ignoring interactions with Fourier harmonics from quad- 
rants other than the first one since they only contribute 
subdominant terms in the short-time expansion of the 
hydrodynamic fields. When such terms are taken into 
account it is likely that singularities obtained at leading 
order will be mostly advected by a modified velocity field 
which carries the singularities slightly out of the y-planc 
without changing their nature, as happens in the work 
of Tanveer and Speziale [18]. Of course, positivity of the 
Fourier coefficients will be lost but not necessarily their 
scaling properties. Observe also that the y-plane being 
a plane of symmetry, this picture implies that there are 
several pieces of the singular manifold very close to the 
y-plane. In Fourier space they produce a kind of inter- 
ference pattern which at first has very long wavelength 
(in k). This wavelength becomes shorter and shorter as 
time advances and the singular manifold moves further 
away from the y-plane. 18 



IV. THE GEOMETRY OF THE 
PSEUDO-HYDRODYNAMIC FLOW 



may rewrite (13)-(14) as a Taylor series in two variables 

oo oo 

^(!/)-2»+!/2=EE^^) Ci 1 C2 2 , (39) 

fei=0fe 2 =0 

d = e^, C2 = e^ 2 . (40) 

If we now hold y 2 (and thus £2) fixed and sum over fc 2 , 
we obtain a Taylor series in £1 such that all its coeffi- 
cients (except possibly the first one) arc positive. By 
Vivanti's theorem [33], if such a series has a finite ra- 
dius of convergence (as is the case here because of the 
aforementioned exponential decrease), the singularity in 
the complex (i-plane nearest to the origin is on the pos- 
itive real axis at a location 2/1 — 2/1(2/2), which depends 
on j/2- The function 2/1 (2/2 ) defines an object which we 
here call the singular manifold and is the edge of the 
analyticity domain 2/1 < 2/1(2/2)- 19 A standard theorem 
about multi-dimensional Taylor series states that their 
domain of convergence is logarithmically convex (see, 
e.g., Ref. [34]). In our case this just means that the 
analyticity domain in the 2/-plane is convex. As shown in 
MBF using slightly different notation, the singular man- 
ifold can be constructed either as the envelope of the 
family of straight lines y\ cos 8 + 2/2 sin 8 = 5(6) (where 
the decrement 8 has been defined in Section III) or as the 
envelope of analyticity disks. 

To numerically construct the pscudo-hydrodynamic so- 
lution in the y-pl&ne from the Fourier data we use (14) 
for the stream function, (18) for the velocity and (20) for 
the vorticity. Although our Fourier data typically have 
35 decimal digits, it suffices to truncate them to 16 dig- 
its to obtain the various relevant fields in j/-space with a 
good accuracy. 



In two-dimensional simulations of hydrodynamics, con- 
siderable insight is usually obtained by looking at flow 
features in the physical space. This is much simpler in 
two dimensions than in three, provided that the relevant 
features are in the real R 2 space. Here the most im- 
portant features are in the complex C 2 space, which is 
equivalent to having four real dimensions. Fortunately, 
as explained in Section II, we can make use of only 
two real dimensions by working in the y-pl&ne above 
(z\, Z2) — (it, 0) which extends in the (pure) imagi- 
nary directions. As already briefly mentioned in Sec- 
tion 4 of MBF, the positivity of all the Fourier coeffi- 
cients F(ki, fc 2 ) (except F(l, 0)) and the exponential de- 
crease with the wavenumber, imply that the solution to 
the (short-time asymptotic) Euler equation has a line of 
singularities S in the (2/1, j/2)-plane. Indeed, since only 
harmonics with non- negative k\ and k 2 are present, we 



A. Presentation of the t/-plane results 

We begin with global topological features and then 
turn to a a more local and more quantitative description. 
Fig. 14 gives a global view of the flow in the y-pl&ne. 20 
The outer edge of the flow region, which passes very 
close to the origin is the singular manifold. At large dis- 
tances on the upper left and the lower right, respectively, 
the singular manifold has logarithmic branches. Close 
to the singular manifold, the streamlines follow it until 
they make a U-turn and eventually plunge into the third 
quadrant (j/i < 0, 2/2 < 0) where they become straight 
with slope 1/2 at large distances. An important feature 
is the Xj-turn separatrix, above which stream lines make 
U-turns which become increasingly sharp when moving 



Somewhat similar interference patterns are obtained when the 
short-time asymptotics is extended to the Navier— Stokes equa- 
tion (with viscosity scaling as 1/t). 



More correctly, the singular manifold is a (perhaps analytic) 
manifold in C 2 whose intersection with the y-plane is designated 
here by the same name. 

When magnifying this figure, ADOBE READER® 7 or higher 
is recommended. 
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FIG. 14: Global geometry of the flow in the y-plane. Stream- 
lines (solid lines) and iso-vorticity lines (thin-dotted lines) are 
shown. Thick-solid-crenated line: singular manifold; thick- 
solid line: U-turn separatrix (ip « 0.5); thick-dashed line: 
vorticity separatrix (to = and ip = In 2). The ticks on 
the two axes correspond to coordinate 0.25. Inset: Contours 
of absolute value of the cotangent of the angle between the 
streamlines and the iso-vorticity lines as a measure of deple- 
tion of nonlinearity. 



to the lower right, and below which there are no U- 
turns. Vorticity contours starting close to the singular 
manifold far on the upper left get pressed increasingly 
close into the singular manifold when moving to the lower 
right. The vorticity separatrix divides negative vorticity 
(above) and positive vorticity (below) . It approaches the 
singular manifold in the lower right but not as fast as the 
U-turn separatrix. In view of the Jacobian formulation of 
the Euler equation, the vorticity separatrix is clearly also 
a streamline. Hence, the strong depletion of nonlinear- 
ity evidenced by accumulation of contour lines near this 
separatrix on the inset of Fig. 14. The depletion is here 
measured by plotting the absolute value of the cotangent 
of the angle between Vi/> and Vw. 

Figs. 15 shows the stream function and the vortic- 
ity with more details in a region of particular interest. 
Increasingly sharp U-turns of the stream lines are seen 
when moving to the lower right into the narrowing chan- 
nel separating the singular manifold from the U-turn sep- 
aratrix. It is seen that the vorticity becomes very large 
and negative near the singular manifold, while the stream 
function remains finite with a value around ip = 0.5, 
the same as on the U-turn separatrix. Thus, the sin- 
gular manifold which is simultaneously a limiting case of 
a streamline and of a vorticity contour, displays strong 
depletion of nonlinearity as seen on the inset of Fig. 14. 
We also looked at the velocity field (not shown); close 
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FIG. 15: Enlargements of Fig. 14 around the point (j/i, j/2) = 
(0.8, —1.0) showing the streamlines (upper figure) and the 
vorticity contours (lower figure). Only negative vorticity con- 
tours are shown. 



to the singular manifold the velocity is parallel to this 
manifold and decreases in modulus when moving down 
and to the right. Note that, contrary to the vorticity, 
the velocity does not grow explosively when approach- 
ing the singular manifold; there is no numerical evidence 
against the plausible assumption that the velocity has a 
finite limit on the singular manifold which is tangent to 
this manifold. Similarly, the pressure (not shown) also 
appears to have a finite limit on the singular manifold. 

For a better quantitative grasp we show in Fig. 16 a 
one-dimensional cut of the two-dimensional fields along 
the normal to the singular manifold passing through the 
origin, shown as a dashed-dotted line on Fig. 14. It is seen 
that the stream function takes the finite value VWg ~ 0.5 
on the singular manifold and the same value at the U- 
turn separatrix and that the vorticity follows approxi- 
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FIG. 16: Upper figure: stream function t/>, velocity compo- 
nents (suitably rescaled) v„ and vt normal and parallel to the 
singular manifold along the line of cut normal to the singu- 
lar manifold passing through the origin, shown as a dashed- 
dotted line on Fig. 14. U and S identify the places where the 
cut intersects the U-turn separatrix and the vorticity separa- 
trix. Lower figure: negative vorticity along the same cut in 
log-log coordinates. 

mately a power law — u oc s _/3 with 0.7 < [3 < 0.9, where 
s is the distance to the singular manifold. The scaling is 
however rather poor; it gets even worse when repeating 
the same analysis along the other dashed line normal to 
the singular manifold, shown on Fig. 14. It is also seen 
that at the singular manifold the normal velocity van- 
ishes linearly. We now turn to comments and theoretical 
explanation of most of these features 



B. Bridging fc-space and y-space results 

As we shall now see, it is quite obvious to relate the 
leading-order asymptotics (23) of the Fourier coefficients 
at large k and the leading-order behavior near the singu- 
lar manifold in y-space. To explain the poor scaling ob- 
served for the vorticity in y-space, we need to take into ac- 
count subleading corrections as we shall also discuss. The 
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FIG. 17: Construction of the singular manifold from the log- 
arithmic decrement 8(0). 

"Fourier-Laplace" representations (14), (18) and (20) for 
the (pseudo-hydrodynamic) stream function, the veloc- 
ity and the vorticity, connect k- and y-space functions. 
Consider, for example, the vorticity; using (23) and polar 
coordinates k = k(cos9, sin6) we can rewrite it as 

u(y) = -^C(6)k- a+2 c- kh{e >y\ (41) 

k 

h(0;y) = S(9)- yi cos9 -y 2 sine . (42) 

The convergence properties at high wavenumbers of this 
sum will depend crucially on the sign of the decrement 
h(6;y). If 

min h(6; y) > , (43) 



all the exponentials are decaying and the sum will be fi- 
nite. If the minimum is negative, the sum is divergent. 
In the borderline case of a vanishing minimum, the alge- 
braic prefactors will determine convergence. If 6(9) is a 
smooth function of 9, as our numerical results suggest, 
the minimum corresponds to a vanishing derivative with 
respect to 9. Hence, the borderline case is characterized 
by the following two equations: 

5(6) -j/i cos 6 -y 2 sin 6 = 0, (44) 
<5'(6>) + yisin6l-y2Cos(9 = , (45) 

where S'(9) is the derivative of 6(6). Those points 
y*(9) = y+ 2 ) which satisfy (44)- (45) are on the sin- 
gular manifold. Conversely, 6(6) is the distance from the 
origin to the tangent at the singular manifold which has 
the slope 6 — ir/2 (see Fig. 17). It follows that the singu- 
lar manifold is the envelope of such lines. In MBF this 
result was derived by Poincare's pinching argument. 

From (44)-(45) it is easily shown that the near-edge 
behavior of 5(9) given by (37)-(38) implies logarithmic 
branches for the singular manifold: y 2 — (1/2) ln(— y\) 
for large negative y\ and y\ ~ ln(— y 2 ) for large nega- 
tive Z/2- 

We observe that s = ming h(6;y) is the shortest Eu- 
clidean distance of y to the singular manifold, with a plus 
sign when y is below the singular manifold and a minus 
sign when it is above. Let us assume that y is below 
or on the singular manifold and let us denote by 6+(y) 



15 



the value of 9 where the minimum is achieved. Near this 
minimum we can Taylor-expand the decrement 



h(0-, y) = s + \k (o - e*f + o{[e- e*f) 



(46) 



where h" = d 2 h(8+, y)/d9 2 . The convergence of the sum 
(41) depends only on the high-A: behavior, where we can, 
to leading order, replace the sum by an integral over 
kdkd9, to obtain a "continuous approximation" 



^cont 



/>27T />00 

(y) = - d9 dkC(9)k~ 
Jo Jo 



a+3 c -kh(6;y) 



■ (47) 



When s vanishes or is small and positive, we can evaluate 
the angular integral in (47) by steepest descent 



fhr f°° 

Wcontfa*) - T77 / dkk ~ 

V n * Jo 



a+5/2 -ks 



(48) 



On the singular manifold, s = and it is clear that the 
integral over k is ultraviolet-divergent as soon as a > 3/2. 
All the values of the prefactor exponent a considered 
in this paper are at least 5/2 and thus give an infinite 
vorticity at the singular manifold. The same analysis 
applied to the stream function and to the velocity gives 
ultraviolet-convergent integrals. For small positive s, we 
obtain from (48) 



Section IV A. Fig. 18 shows — w synt h (y) in log-log coor- 
dinates for three values of n. The lowest one n = 1000 is 
comparable to what is used in the actual Euler calcula- 
tion: the scaling is very poor. Only when we increase the 
resolution more than tenfold to n = 11, 000 do we begin 
to see a clean s -5 / 6 scaling. 
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FIG. 18: Same as lower part of Fig. 16 but for the synthetic 
data with various values of the maximum wavenumber n. In- 
set: corresponding logarithmic local slopes. The predicted 
scaling exponent is —5/6. 



uj cont (y*) ~ -0(0+) J — 1X7/2-0)*-" , (49) 



2~ Q 



(50) 



where T(-) denotes the Gamma function. 

For SOC initial conditions, a s=s 8/3 and thus the vor- 
ticity diverges to leading order with a s -5 / 6 law, when 
approaching the singular manifold. The subleading cor- 
rections causing the poor scaling seen in Section IV A are 
of various sorts. First, there are subleading corrections 
to (23) whose simplest manifestation is the discrepancy 
between the local scaling exponent a\ oc (k) and its extrap- 
olated value aoo, as discussed in Section III A. As already 
stated, we do not know the functional form of such cor- 
rections. Second there are subleading corrections coming 
from having approximated the Fourier-Laplace sums by 
integrals. It is easily shown that they contribute O(s ) to 
the vorticity in y-space. Third there are the subleading 
corrections to the continuous approximation (47), which 
may be shown to be 0(s~ 1 ^ 3 ). A simple way to deter- 
mine how much the scaling is degraded by the second 
and third type of corrections is to usesynthetic data for 
which the Fourier coefficients are given exactly by (23). 
It is then easy to change the resolution n x n and to find 
out how large n should be for clean leading-order scaling 
to emerge. We performed such a calculation with C and 
a constant and the values of 5(9) taken from the actual 
Euler SOC data. From the synthetic data, using (20) we 
then calculate a synthetic vorticity w syn th(y) just as in 



C. Theory of j/-plane pseudo- hydrodynamics 

The starting point for theory in the y-plane is of 
course the pseudo-hydrodynamic vorticity equation and 
its boundary conditions far into the third quadrant, de- 
rived in Section II and repeated here for convenience: 



v ■ Vw + uj = 0, 



(51) 



v ~ (-1, -0 , u~ ic^ - 2e 2 ^, Vl , y 2 -> -oo.(52) 

Alternatively we can rewrite (51) as J(ip,u>) — lu or as 
J(ip, In \u>\) = 1. Thus the map from (y\, y 2 ) to (ip, In |w|) 
is area preserving. The vorticity separatrix was defined 
by u — 0; so that that the Jacobian of the stream func- 
tion and of the vorticity is zero along this line. Hence it 
is also a streamline. 21 It is easily shown that the value 
of the streamfunction on this line is In 2. Indeed, as we 
follow the vorticity separatrix far into the third quad- 
rant, we obtain from (52) that y\ ~ 2y 2 + 2 In 2. Since 
ip = (1/2)?/! — y 2 (up to exponentially small terms), we 
obtain the result claimed. We have also checked numeri- 
cally that the value of the stream function on the vorticity 
separatrix is In 2 to at least three decimal places. 



21 By changing u) into in (51), we can show similarly that the 
singular manifold, at which = 0, is also a streamline. 
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Depletion of nonlinearity near the singular manifold 
prevents us from using the dynamical equation (51) to 
derive the scaling exponent of the singularities by, e.g., 
balancing to leading order the two terms in (51). Nev- 
ertheless, such balancing gives some useful information, 
such as the vanishing of the normal component v n (s) of 
the velocity near the singular manifold (for s — > 0) and 
the independence on position of the exponent (3 charac- 
terizing the divergence of the vorticity. Since (3 and the 
prefactor exponent a are related by a + j3 = 7/2, this will 
establish the independence of a on 9, which was rather 
strongly supported by the numerical results reported in 
Section III. We now derive these results. In what follows, 
points on the singular manifold are parameterized by 
the angle 9 between the yi-axis and the outgoing normal. 

To show the vanishing of v n (s) for s — ► 0, it is conve- 
nient to use as local coordinates near the singular mani- 
fold the angle 9 and the distance s. We denote by v n (s, 9) 
and Vt(s, 9) the components of the velocity along the in- 
ward normal and along the tangent in the direction of 
increasing 9. For small s, to leading order, (51) becomes 




v n (s, 9)d s + v t (s, 9)— \— d e 



u(s,6) ~ -uj(s,6) , (53) 



FIG. 19: A typical streamline passing near the singular man- 
ifold and performing a U-turn. Although this figure uses the 
actual data rather then being a sketch, the distance between 
the streamline and the singular manifold has been somewhat 
increased for legibility. 



the smallest velocities and thus the leading-order con- 
tribution to the integral J^ 1 d£/\v(£)\ are expected to 
come from the immediate neighborhood of the U-turn. 
We have checked this conjecture numerically by calculat- 
ing the vorticity and the modulus of the velocity along 
a streamline chosen to have a rather sharp but well- 
resolved U turn (see Fig. 20). 



where R(9) is the radius of curvature of the singular man- 
ifold (the arclength is given by R(9)d9). We now assume 
that u) oc s _/3 (with (3 > 0). If v n (0, 6) did not vanish, 
the first term on the l.h.s. of (53) would be proportional 
to s - ^ -1 , which for small s could not be balanced by any 
of the other terms of the equation. The argument actu- 
ally implies the stronger result that the normal velocity 
v n cannot vanish more slowly than s 1 . The technique 
of Section IV B on bridging fc-space and y-space results 
can be used to show that dv n /ds remains finite at the 
singular manifold, although it has the same dimension as 
the vorticity which becomes infinite. 22 Thus v n actually 
vanishes linearly with s. 

For the independence of (3 on 9, we integrate (51) along 
a typical streamline passing near the singular manifold 
between two points Mq (far from the singular manifold) 
and Mi (within a small distance si), so that there is a 
U-turn in between which is assumed not to be close to 
either M or Mi (see Fig. 19). We obtain 
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FIG. 20: Modulus of the velocity and vorticity (changed sign) 
vs arclength I along the streamline shown in Fig. 19 which has 
ip ~ 0.47. The point Mo is taken near I — 2.5; the U-turn, 
Mi and M2 are near I = 3, I = 3.2, i — 3.5, respectively. 



a; (Mi) 
w(M ) 




The streamline is here parameterized by the arclength I 
measured from an arbitrary reference point and growing 
when moving into the upper far left, opposite to the di- 
rection of the velocity. When moving from Mq to Mi 



In the proof one uses the vanishing of the normal component 
of the velocity at the singular manifold, a consequence of the 
singular manifold being a streamline. 



We assume now that the vorticity near the singular 
manifold is given to the leading order by u> — c(6)s~^^\ 
where we temporarily leave the possibility that the ex- 
ponent (3 depends on the parameter 9 associated to the 
nearest point on the singular manifold. We take a second 
point Mi on the same streamline but further away from 
the U-turn. We then have (to leading order) 

lo(Mi) ~ c(0i)*r Wl) , (55) 

w(M 2 ) ~ c(0 2 )s2 W2) . (56) 

where (si, 0\) and (S2, 6*2) are the local coordinates for 
Mi and M 2 . From (54), applied successively to Mi and 
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M 2 we find that 

w(M 2 ) 



(M; 



exp 




(57) - 



Between Mi and M 2 the streamline is close to the singu- 
lar manifold and the velocity is dominated by its tangen- 
tial component; hence we can replace the r.h.s. of (57) 
by an integral along the singular manifold and obtain to 
leading order 



uj{M 2 ) 
w(Mi) 



K\2 = exp 




(58) 



which depends neither on si nor on s 2 . We now ob- 
serve that the solenoidal character of the velocity implies 
(again to leading order) 



siv t (6i) ^ s 2 v t {9 2 ) 
It follows from (58) and (59) 

-13(02) 



oj(M 2 ) ~ c 2 



SlVl 

v 2 



K 12 c 1 s l 



(59) 



(60) 



Comparison of the middle and the rightmost members 
gives 



(3{d 1 )=(3{6 2 ), c 2 = K 12Cl 



(61) 



This establishes the independence of the vorticity scaling 
exponent on 9. 



V. A PASSIVE SCALAR MODEL 

As shown in Ref. [35] simple advection of a passive 
scalar by a prescribed velocity field with just a few 
Fourier harmonics can easily lead to singularities because 
fluid particles may come from or go to (complex) infin- 
ity in a finite time. In the present context of short-time 
asymptotics, the equivalent of a passive scalar model is 
to treat the (pseudo-hydrodynamic) vorticity u> in (21) as 
a passive scalar advected by a prescribed velocity. The 
simplest prescribed velocity we can take is 



vp(y) 



e 2y \ ^c yi 



obtained from the stream function 

1 1 1 



2 V1 



2/2 



2 



2 6 



2y 2 



(62) 



(63) 



This velocity field includes the drift (— 1, —1/2) result- 
ing from the shifts of the original coordinates by terms 
proportional to In t and the contributions from the basic 
modes. For our passive scalar model we use the vortic- 
ity equation (21) with the inclusion of an inhomogeneous 




FIG. 21: Streamlines in the y plane for the passive scalar 
model given by (63), which has a hyperbolic stagnation point 
at the origin. Thick line with arrows pointing to the origin: 
stable manifold of the associated dynamical system (65) which 
is also the singular manifold for the vorticity. Thick line with 
arrows pointing away from the origin: unstable manifold. 



term whose precise form does not matter (as long as it 
does not have itself any singularity): the singularities of 
the passive vorticity stem solely from advection. Specif- 
ically, the passive scalar model is defined by 



V P -X7LU + LU=-C yi+2y2 



(64) 



where the r.h.s. is taken to be the interaction term of the 
two basic modes. It is easy to write down Fourier-space 
recursion relations for this model and to show that all 
the Fourier coefficients are positive. Numerical solution 
of the recursion relations gives the usual type of scaling 
with a very clean prefactor exponent a = 5/2. 23 In the 
y-space this implies a blow up of the vorticity w a r 1 
as function of the distance s to the singular manifold. 

Actually all these results can be derived in a rather 
straightforward manner by working in the y-space, as we 
now explain. Eq. (64) can be integrated along the charac- 
teristics. For this we consider the conservative dynamical 
system of fluid particle trajectories in the velocity field 



d_ 



V = v P (y) . 



(65) 



The integral lines are the lines ipp = const., which are 
shown in Fig. 21. At the origin there is a hyperbolic 
stagnation point near which we have ipp = — (l/4)yf + 
0(\y\ 3 ). The associated unstable manifold is simply 



23 With the already cited asymptotic interpolation method [31] the 
exponent a is found to differ from 5/2 by less than 10 -11 . 
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yi = 2j/ 2 , while the stable manifold is the other solution 
to V'p = 0. 

This hyperbolic stagnation point at the origin com- 
pletely determines the scaling of the vorticity. Indeed, in 
Section IV C we derived from the vorticity equation (21) 
an expression (54) which shows that large vorticities stem 
from low- velocity regions. This derivation, which did not 
make use of the fact that the vorticity is the curl of the ve- 
locity, remains valid for the passive scalar model (except 
for minor changes due to the presence of an inhomogc- 
neous term). In the fully nonlinear case, the low veloci- 
ties were due to the increasingly sharp U-turns described 
in Section IV A. In the present much simpler case, they 
are just due to the passage near the hyperbolic stagna- 
tion point. As we follow a streamline upstream (i.e. to 
increasing arclengths I in the notation of Section IV C) 
we come closer and closer to the stable manifold. The 
latter thus plays the role of the singular manifold. By 
the same argument as used in Section IV C the scaling 
of the vorticity near the stable/singular manifold is the 
same everywhere. It suffices to determine it locally near 
the stagnation point. One way is to parameterize the 
streamline by the distance s to the stable manifold (more 
precisely to its tangent at the origin, since we are doing 
a local analysis). By (54), the growth of the vorticity is 
controlled by 

exp {/R^} =cxp {/^)}' (66) 

where v n (s) is the velocity component along the (inward) 
normal n = (— 1/V5, — 2/\/5) to the stable manifold at 
the origin. Near a hyperbolic stagnation point we have 
v n (s) = As + 0(s 2 ), where A is the positive eigenvalue of 
the velocity gradient diVj at the stagnation point. Here, 
it is elementary to show that A — 1. Using this in (66) 
we find that the vorticity ui oc s" 1 , as claimed. This 
argument is easily adapted to the passive scalar mod- 
els for other two-mode initial conditions, such as those 
discussed in Section IIIC. The same s _1 behavior is al- 
ways obtained, which is thus universal, contrary to what 
happens in the full nonlinear case. 

Although the passive scalar model does not predict the 
exact and non-universal character of the vorticity blow 
up for the full nonlinear problem, the singular manifold 
is given quite accurately by the passive scalar model. 
In particular it is immediately checked that both have 
the same logarithmic branches (at least to leading or- 
der). Furthermore in the passive scalar model the sta- 
ble/singular manifold goes exactly through the origin 
while in the full nonlinear problem it passes within a dis- 
tance 8 « 0.0065, as shown in MBF. It may be that such 
agreements are due to the presence of very strong deple- 
tion of nonlincarity in the full problem, thereby making 
a simple linear advection model quite relevant. 

Actually, the passive scalar model can be systemati- 
cally improved by enriching the prescribed velocity field 
through addition of higher-order modes. A simple way 



to do this is to take all the Fourier modes such that 
ki + < n + 1. For SOC, we have studied these 

"enriched" passive scalar models for various values of n. 
They all possess a hyperbolic stagnation point. The asso- 
ciated positive eigenvalue A„ becomes larger than unity 
when n > 1. The first few values for the correspond- 
ing prefactor exponent a n = 7/2 — 1/A„ arc: a = 2.5, 
at\ s=s 2.594, a 2 ~ 2.613. For larger values of n the growth 
is very slow; for example a 2 o ~ 2.618. We also observed 
that as n grows, the stagnation point moves to the right 
and down and the angle between its stable and unsta- 
ble manifolds decreases. It is likely that for n — > oo, 
the stagnation point is pushed to infinity in such a way 
that its stable and unstable manifold tend to the singular 
manifold and to the U-turn separatrix for the nonlinear 
problem, while a n — > a, but the convergence may be 
slow. 



VI. CONCLUSION 

The present paper, like Refs. [4] and MBF, is mainly 
concerned with the short-time asymptotics of the 2D 
Euler equation in situations where complex-space sin- 
gularities are born at infinity at time t = 0+. Let us 
first summarize the main findings of this work, which 
uses a mixture of ultra-high precision computations (with 
up to 100-digit accuracy) and of theory Our work is 
specifically concerned with initial conditions in the form 
of a trigonometric polynomial; it is shown in the Ap- 
pendix that this problem can generically be reduced to 
one with only two modes. A very detailed description of 
the complex singularities is given. For all cases studied, 
the Fourier coefficients except one are found to be non- 
negative (this was already reported for SOC in MBF). In 
any direction of rational slope tan 9 not too close to the 
edges of the Fourier domain, the coefficients of the stream 
function converge very quickly with increasing wavenum- 
bers k to the form C(9)k~ a e~ kS( - 9 \ The prefactor expo- 
nent a, determined with better than one percent accu- 
racy, is independent of 9 but is not universal: when the 
initial modes are orthogonal, it is indistinguishable from 
8/3 w 2.66, whereas with a 45 degree angle between the 
initial modes it takes the value 2.54. We cannot rule out 
that a depends also on the moduli of the initial modes 
but we have no evidence that it does. 

It is shown that the singularity problem can 
be reformulated as an ordinary steady-state 
(pseudo)hydrodynamic problem in a suitable y-plane 
corresponding to pure imaginary coordinates. The 
complex singularities are in this y-plane on a smooth 
(possibly analytic) curve extending to infinity with 
logarithmic branches. The vorticity diverges as s - ' 3 , 
where s is the distance to the singular manifold and 
a + [3 = 7/2. We give a full description of the geometry 
of streamlines and vorticity contours in the y-plane 
(Fig. 14). Increasingly sharp U-turns of the streamlines 
near the lower logarithmic branch of the singular 
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manifold give rise to the vorticity scaling. Very strong 
depletion of nonlincarity near the singular manifold 
prevents application of dominant balance to determine 
the scaling exponent of singularities and is likely to 
be the reason for the very unusual non-universality of 
the singularities. Finally it is shown that the scaling 
behavior of the prefactor persists in time significantly 
beyond the validity of the short-time asymptotics, at 
least as intermediate asymptotics. However, we do not 
know if the non-universality of the singularities found 
in the short-time regime carries over to the full Eulcr 
equation. 

The main theoretical shortcomings of this work are our 
inability so far to prove the positivity of Fourier coeffi- 
cients and to derive the prefactor exponent a (or the 
vorticity divergence exponent (3) from the initial condi- 
tions (we also failed to identify the nature of subleading 
corrections to (23)). We have nevertheless gained some 
qualitative understanding with the passive scalar model 
of Section V that ignores the back reaction of the vortic- 
ity on the velocity but which sheds interesting light on 
the mechanism for producing singularities. In this toy 
model, scaling is controlled by a stagnation point of the 
velocity field, whereas in the full nonlinear problem the 
stagnation point is rejected to infinity. 

We have described our findings in some detail, hoping 
that colleagues will be able to help us with the missing 
theory. 

In principle the methods used for the 2D short-time 
Euler problem can be extended to various other short- 
time problems. One instance is the short-time regime 
for the 2D ideal incompressible MHD equations. A pre- 
liminary study for this case indicates that the positivity 
result does not survive: Fourier amplitudes display oscil- 
lations revealing a richer geometry of the singular mani- 
fold which can no more be captured in terms of just the 
imaginary coordinates y. In mathematical terms, one 
has to study the amoeba and coamoeba of the singular 
manifold. 24 Oscillations can be handled by techniques 
similar to those discussed here, as has already been done 
in Rcf. [19]. 

Another natural extension of our study is to the 3D Eu- 
ler equations which also have a short-time regime. This 
is rather straightforward. A direct extension of the al- 
gorithm used in two dimension requires CPU resources 
(time complexity) proportional to fc^ ax instead of fc^ ax - 
This becomes prohibitively large when fc max exceeds a 



In d-dimcnsional algebraic geometry one deals with an algebraic 
manifold in complex coordinates Cl > • • ■ Cd an d the amoeba is de- 
fined as the image of the manifold under the map Cl 1 — ^ J/i = 
In . . . ,fd I— > y d = ln|fd|- The coamoeba is similarly de- 
fined in terms of the argument functions of the f 's. The com- 
plex exponentials e" 1 Z1 , . . . , e" 1 Zd , play here the role of the Q's. 
The name amoeba has been proposed by Gelfand, Kapranov and 
Zelevinsky [36] because amoebae sometimes have pseudopods re- 
sembling those of microscopic protozoa. Coamoebae have been 
introduced by Tsikh and Passare (private communication). 



few hundred. In principle, the time complexity can be 
reduced to fc^ ax (with logarithmic corrections) by using 
FFT's and the recent technique of "relaxed multiplica- 
tions" [37] . However in the calculations reported in MBF 
and the present paper the magnitude of Fourier coeffi- 
cients can vary by several hundred orders of magnitude; 
this requires special precautions when applying FFT's 
unless one is prepared to use several hundred digits. 

We remind the reader that our long term goal is to 
find out about blow up in three dimensions (3D). We 
hope this will not take another 250 years. Progress may 
however be painfully slow if, as we expect, numerical ex- 
perimentation is to play an important part. Indeed, the 
amazingly fast growth of computer power observed over 
the last 50 years, becomes much less spectacular when 
translated in terms of resolution achievable in 3D simu- 
lations. 25 . As more powerful computers become available 
for investigation of 3D blow up, it would not be advis- 
able to use the new resources exclusively for increasing 
the spatial resolution. Experience on the advantage of 
ultra-high precision for singularity studies from the work 
of Krasny [21], of Shelley [20], of Caflisch [19] and also 
from our own work suggest that it is not safe to use less 
than 30-35 digits. Using flows with symmetries such as 
the Taylor-Green [13, 14] or the Kida-Pelz [15-17] flow to 
boost the resolution introduces a possible element of non- 
genericity, but we can always use such flows to sharpen 
our tools and then, as computers become more powerful, 
turn to flows without symmetry. 

Have the results reported in MBF and the present pa- 
per brought us closer to this Holy Grail of 3D blow-up? 
In a direct way, we cannot infer anything regarding 3D 
real blow-up from a 2D study of complex singularities 
at short times. We have however learned that in this 
rather restricted framework, singularities are located on 
very smooth objects (possibly analytic manifolds); be- 
cause the fastest spatial variation is then in the direction 
perpendicular to the singular manifold, the singularities 
have strongly depleted nonlinearity in a suitable frame. 
We have already good evidence that in 2D this smooth- 
ness property is not limited to the short-time regime [4]. 
In 3D such a property would be both a curse, since dom- 
inant balance cannot be used, and perhaps a blessing, 
since it might well slow down (indefinitely?) the ap- 
proach of singularities to the real domain. 

Acknowledgments 

We are grateful to M. Blank, H. Frisch, J. van der Ho- 
even, D. Mitra, A. Pumir, A. Sobolevskh, A.K. Tsikh, 
P. Zimmcrmann and an anonymous referee for useful 
discussions and comments. Part of this work was done 
while the authors participated in "Frontiers of Non Lin- 



25 At the moment the highest resolution accessible in 15-digit pre- 
cision for three-dimensional flow without any special symmetries 
is 2048 3 [38, 39] 



ear Physics" (Nizhny Novgorod, Russia, July 5-12, 2004) 
and in "Singularities, coherent structures and their role 
in intermittent turbulence" (Warwick, UK, September 
9-17, 2005). TM and JB were supported by the Grant- 
in- Aid for the 21st Century COE "Center for Diversity 
and Universality in Physics" from the Japanese Ministry 
of Education. TM was supported by the Japanese Min- 
istry of Education Grant-in- Aid for Young Scientists [(B) , 
15740237, 2003] and by the French Ministry of Educa- 
tion. UF was supported by the Grant-in-Aid for the 
21st Century COE "Center of Excellence for Research 
and Education on Complex Functional Mechanical Sys- 
tems" from the Japanese Ministry of Education. WP 
had partial support from the European network EU RTN 
no. HPRN-CT-2002-00282 'Hyke. JB, UF and WP had 
partial support from the Federation de Recherche Wolf- 
gang Docblin (FR 2800 CNRS). Part of the computa- 
tional resources were provided by the Yukawa Institute 
for Theoretical Physics (Kyoto) and by the Mesocentre 
SIGAMM (Nice). 



APPENDIX A: REDUCTION OF MULTIMODE 
INITIAL CONDITIONS 

Here we shall show that the short-time asymptotics of 
two-dimensional Euler flows with generic initial condi- 
tions of trigonometric polynomial type can be reduced to 
the study of two-mode initial conditions. In Section II 
we have seen that with two initial modes p and q the be- 
havior of the stream function ty(z,t) for large imaginary 
arguments |yi|, 1 2/2 1 can be described by the similarity 
ansatz (6)- (7). It relies on the fact that when y is such 
that P2/P1 < 2/2/2/1 < I2/I1 the leading-order factors ac- 
companying each factor t in the time- Taylor expansion 
(4) are either e~ Jpz or e~* qz . In the limit \y\ — > 00, 
t^Owe can make such terms finite by shifting simulta- 
neously p ■ y and q • y by Int. 

The similarity ansatz as explained above is however 
not applicable to the case of more than two initial modes. 
Instead, we have to reduce the multimode initial condi- 
tion to various two-mode problems which can be handled 
in the usual way. Let us illustrate this by looking at a 
simple three-mode initial condition 

*o(as) = foe**" + h 2 e iqx + +h 3 e"- x + c.c. , (A.l) 

in which the vectors p, q and r are listed in angular coun- 
terclockwise order. As in Section II, to avoid pathologies, 
we assume that the vectors p, q, x are not parallel and not 
of the same length. In the Taylor expansion (4) each fac- 
tor t will now be accompanied by a factor e~* pz , e~ iq z 
or e~ ir ' z . Note that, in the limit \y\ — > 00, t — > we can- 
not simultaneously make the terms te py , te q v and te x y 
remain finite. Indeed, we cannot translate y in such a 
way that all three scalar products p ■ y , q • y and t • y are 
shifted by Int. If P2/P1 < 2/2/2/1 < Wli the factors te Vy 
and te qy will dominate, while if q2/qi < 2/2/2/1 < ^2/^1 
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k. 




FIG. 22: Construction of relevant Fourier modes in various 
angular sectors. Black circles: initial modes; dash-dotted 
lines: boundary of the Newton polytope of the initial modes; 
white circles: higher-order modes associated to nth genera- 
tion. Thick lines show the edges of the various angular sec- 
tors. 



the factors te qy and te ry will dominate. In each case 
the three-mode initial condition (A.l) is reduced to a 
two-mode problem involving either p and q or q and t. 

Let us now turn to the general multimode case with 
an initial stream function of the form 

(fc 1 ,fc 2 )esupp_F(°) 

(A.2) 

Here, we assume Hermitian symmetry 26 F^*{k\,k 2 ) = 
F(°\—ki,—k 2 ) and we take the sum over all wavevec- 
tors for which the Fourier coefficients F^°\ki, k 2 ) do not 
vanish, called the support of F^ and denoted suppi^ ). 
The fact that supp F^ is a finite set plays a crucial part 
in our analysis. Let us suppose for simplicity that all 
initial modes have different lengths, that is K^i,^)! 7^ 
\(k'{,k%)\ for all pairs {k[,k' 2 ) and {k'{,k' 2 ') € suppF^. 

As we have seen before, it is necessary to distinguish 
between different directions in the y-space when taking 
the limits \y\ — > 00, t — > 0. Therefore we let \yi\, \y 2 \ — ► 
00 while keeping the ratio 2/2/2/1 fixed. By Hermitian 
symmetry, it is enough to consider only the case 2/1 — * 
+00. Assuming as in Section II that 



In fact, this condition is not essential. We could as well consider 
more general sets supp F^ ' . 
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*(z,t) = Y,* n (z)t n , (A.3) 

n=0 

and denoting by (k) the Fourier coefficients of \I/ n we 
obtain easily from the Euler equation the following recur- 
sion relations for the Fourier coefficients of the (n + l)th 
"generation" : 

n+i l K l m+p=„fc'+fe"=fe (A.4) 

(fe'Afe")|fe"| 2 -F( m )(fc;,^)FW(^',^') , 

which allows us to compute F^ n+1 \k) in terms of the 
previous generations F^ m \k) and p( p \k) with m,p > 
and m + p = n. From (A.4) follows immediately that 
has finitely many non- vanishing modes 
We now identify the modes in the nth generation which 
give the leading-order contributions to ^ n {z) for a fixed 
ratio 2/2/2/1- For this we use the notion of Newton poly- 
tope of suppF' ). It is defined as the convex hull (in the 
usual sense) of the set suppF^ \ as for example repre- 
sented (taking into account the Hermitian symmetry) on 
Fig. 22 for a typical initial condition. We shall call rel- 
evant those initial modes lying on the boundary of the 
Newton polytope of supp F^ (black circles indicated by 
arrows on Fig. 22). The relevant modes divide the k- 
space into angular sectors, analogously to the three-mode 
case presented above. 



Let now k' and k" be two relevant modes defining an 
angular sector such that k' A k" > 0. The vectors k' 
and k" define an angular sector in the y-space such that 
k'2/k'i < 2/2/2/1 < k'2/k" . For these directions the leading- 
order terms in ^ n (z) are proportional to e^" k +n k s> ' v 
with n', n" > 1 and n'+n" = n+1; the terms correspond- 
ing to other Fourier modes are subdominant. Clearly, we 
would have obtained the same dominant terms if we had 
started from just the two initial modes k' and k" . We 
can now apply the similarity ansatz, shifting k' ■ y and 
k" ■ y by Int. Let us remark that while it is possible to 
eliminate the time variable by a global similarity ansatz 
for two-mode asymptotics, it is in general impossible to 
do so for more than two modes. 

In the exceptional cases where two or more relevant 
modes have the same length one must take into account 
the fact that the Fourier coefficient of their sum vanishes. 
The description of the leading-order contributions in the 
nth generation becomes slightly more involved than in 
the generic case, but the leading-order behavior in a fixed 
direction 2/2/2/1 is still dominated by two-mode asymp- 
totics. 

Summarizing, we have shown that the study of the 
short-time asymptotics of Euler flows with multimode 
initial conditions can be reduced to the analysis of various 
two-mode asymptotics in angular sectors defined by a 
suitable set of relevant initial modes, namely those on the 
boundary of the Newton polytope of the initial modes. 
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